{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using VMLS\n",
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 7\n",
    "# Matrix examples\n",
    "### 7.1 Geometric transformations\n",
    "\n",
    "Let’s create a rotation matrix, and use it to rotate a set of points $π/3$ radians ($60^{\\circ}$). The result is in Figure 7.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Array{Float64,2}:\n",
       " 0.5       -0.866025\n",
       " 0.866025   0.5     "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Rot(theta) = [cos(theta) -sin(theta); sin(theta) cos(theta)];\n",
    "R = Rot(pi/3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"500\" height=\"500\" viewBox=\"0 0 2000 2000\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip1500\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2000\" height=\"2000\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip1500)\" points=\"\n",
       "0,2000 2000,2000 2000,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip1501\">\n",
       "    <rect x=\"400\" y=\"200\" width=\"1401\" height=\"1401\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip1500)\" points=\"\n",
       "153.898,1887.47 1952.76,1887.47 1952.76,47.2441 153.898,47.2441 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip1502\">\n",
       "    <rect x=\"153\" y=\"47\" width=\"1800\" height=\"1841\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip1502)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  235.664,1887.47 235.664,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1502)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  644.496,1887.47 644.496,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1502)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1053.33,1887.47 1053.33,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1502)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1462.16,1887.47 1462.16,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1502)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1870.99,1887.47 1870.99,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1502)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  153.898,1803.83 1952.76,1803.83 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1502)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  153.898,1385.59 1952.76,1385.59 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1502)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  153.898,967.359 1952.76,967.359 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1502)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  153.898,549.125 1952.76,549.125 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1502)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  153.898,130.891 1952.76,130.891 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1500)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,1887.47 1952.76,1887.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1500)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,1887.47 153.898,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1500)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  235.664,1887.47 235.664,1859.87 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1500)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  644.496,1887.47 644.496,1859.87 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1500)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1053.33,1887.47 1053.33,1859.87 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1500)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1462.16,1887.47 1462.16,1859.87 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1500)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1870.99,1887.47 1870.99,1859.87 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1500)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,1803.83 180.881,1803.83 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1500)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,1385.59 180.881,1385.59 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1500)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,967.359 180.881,967.359 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1500)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,549.125 180.881,549.125 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip1500)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  153.898,130.891 180.881,130.891 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip1500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 235.664, 1939.47)\" x=\"235.664\" y=\"1939.47\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 644.496, 1939.47)\" x=\"644.496\" y=\"1939.47\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1053.33, 1939.47)\" x=\"1053.33\" y=\"1939.47\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1462.16, 1939.47)\" x=\"1462.16\" y=\"1939.47\">1.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1870.99, 1939.47)\" x=\"1870.99\" y=\"1939.47\">2.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 133.898, 1821.33)\" x=\"133.898\" y=\"1821.33\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 133.898, 1403.09)\" x=\"133.898\" y=\"1403.09\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 133.898, 984.859)\" x=\"133.898\" y=\"984.859\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 133.898, 566.625)\" x=\"133.898\" y=\"566.625\">1.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip1500)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 133.898, 148.391)\" x=\"133.898\" y=\"148.391\">2.0</text>\n",
       "</g>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1053.33\" cy=\"1803.83\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1053.33\" cy=\"1803.83\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1462.16\" cy=\"1803.83\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1462.16\" cy=\"1803.83\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1870.99\" cy=\"1803.83\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1870.99\" cy=\"1803.83\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1053.33\" cy=\"1594.71\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1053.33\" cy=\"1594.71\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1462.16\" cy=\"1594.71\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1462.16\" cy=\"1594.71\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1053.33\" cy=\"1385.59\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1053.33\" cy=\"1385.59\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"644.496\" cy=\"1079.42\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#e26f46; stroke:none; fill-opacity:1\" cx=\"644.496\" cy=\"1079.42\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"848.911\" cy=\"717.223\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#e26f46; stroke:none; fill-opacity:1\" cx=\"848.911\" cy=\"717.223\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1053.33\" cy=\"355.022\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#e26f46; stroke:none; fill-opacity:1\" cx=\"1053.33\" cy=\"355.022\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"467.466\" cy=\"974.866\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#e26f46; stroke:none; fill-opacity:1\" cx=\"467.466\" cy=\"974.866\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"671.882\" cy=\"612.665\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#e26f46; stroke:none; fill-opacity:1\" cx=\"671.882\" cy=\"612.665\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"290.437\" cy=\"870.308\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip1502)\" style=\"fill:#e26f46; stroke:none; fill-opacity:1\" cx=\"290.437\" cy=\"870.308\" r=\"14\"/>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create a list of 2-D points\n",
    "points = [ [1,0], [1.5,0], [2,0], [1,0.25], [1.5, 0.25], [1,.5] ];\n",
    "# Now rotate them.\n",
    "rpoints = [ R*p for p in points ];\n",
    "# Show the two sets of points.\n",
    "using Plots\n",
    "scatter([c[1] for c in points], [c[2] for c in points])\n",
    "scatter!([c[1] for c in rpoints], [c[2] for c in rpoints])\n",
    "plot!(lims = (-0.1, 2.1), size = (500,500), legend = false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure 7.1** Counterclockwise rotation by 60 degrees applied to six points.\n",
    "\n",
    "### 7.2 Selectors\n",
    "**Reverser matrix.** The reverser matrix can be created from an identity matrix by reversing the order of its rows. The Julia command reverse can be used for this purpose. `(reverse(A,dims=1)` reverses the order of the rows of a matrix; `flipdim(A,dims=2)` reverses the order of the columns.) Multiplying a vector with a reverser matrix is the same as reversing the order of its entries directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "reverser (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reverser(n) = reverse(eye(n),dims=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 0.0  0.0  0.0  0.0  1.0\n",
       " 0.0  0.0  0.0  1.0  0.0\n",
       " 0.0  0.0  1.0  0.0  0.0\n",
       " 0.0  1.0  0.0  0.0  0.0\n",
       " 1.0  0.0  0.0  0.0  0.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = reverser(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Float64,1}:\n",
       " 5.0\n",
       " 4.0\n",
       " 3.0\n",
       " 2.0\n",
       " 1.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = [1., 2., 3., 4., 5.];\n",
    "A*x # Reverse x by multiplying with reverser matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Float64,1}:\n",
       " 5.0\n",
       " 4.0\n",
       " 3.0\n",
       " 2.0\n",
       " 1.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reverse(x) # Reverse x directly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Permutation matrix.** Let’s create a permutation matrix and use it to permute the entries of a vector. In Julia, there is no reason to create a matrix to carry out the permutation, since we can do the same thing directly by passing in the permuted indexes to the vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       "  2.4\n",
       "  0.2\n",
       " -1.7"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [0 0 1; 1 0 0; 0 1 0]\n",
    "x = [0.2, -1.7, 2.4]\n",
    "A*x # Permutes entries of x to [x[3],x[1],x[2]]\n",
    "x[[3,1,2]] # Same thing using permuted indices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.3 Incidence matrix\n",
    "**Incidence matrix of a graph.** We create the incidence matrix of the network shown in Figure [7.3](https://web.stanford.edu/~boyd/vmls/vmls.pdf#figure.7.3) in VMLS. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×5 Array{Int64,2}:\n",
       " -1  -1   0   1   0\n",
       "  1   0  -1   0   0\n",
       "  0   0   1  -1  -1\n",
       "  0   1   0   0   1"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " A = [ -1 -1 0 1 0; 1 0 -1 0 0 ; 0 0 1 -1 -1 ; 0 1 0 0 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Int64,1}:\n",
       "  1\n",
       " -1\n",
       "  1\n",
       "  0\n",
       "  1"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xcirc = [1, -1, 1, 0, 1] # A circulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Array{Int64,1}:\n",
       " 0\n",
       " 0\n",
       " 0\n",
       " 0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A*xcirc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Array{Float64,1}:\n",
       " 1.1102230246251565e-16\n",
       " 0.0                   \n",
       " 0.0                   \n",
       " 0.0                   "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s = [1,0,-1,0]; # A source vector\n",
    "x = [0.6, 0.3, 0.6, -0.1, -0.3]; # A flow vector\n",
    "A*x + s # Total incoming flow at each node"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Dirichlet energy.** On page [135](https://web.stanford.edu/~boyd/vmls/vmls.pdf#section*.172) of VMLS we compute the Dirichlet energy of two potential vectors associated with the graph of Figure [7.2](https://web.stanford.edu/~boyd/vmls/vmls.pdf#figure.7.2) in VMLS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×5 Array{Int64,2}:\n",
       " -1  -1   0   1   0\n",
       "  1   0  -1   0   0\n",
       "  0   0   1  -1  -1\n",
       "  0   1   0   0   1"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [ -1 -1 0 1 0 ; 1 0 -1 0 0 ; 0 0 1 -1 -1; 0 1 0 0 1 ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Array{Int64,1}:\n",
       " 1\n",
       " 2\n",
       " 2\n",
       " 1"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vsmooth = [ 1, 2, 2, 1 ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Array{Int64,1}:\n",
       "  1\n",
       " -1\n",
       "  2\n",
       " -1"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vrough = [ 1, -1, 2, -1 ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.9999999999999996, 27.0)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(A'*vsmooth)^2, norm(A'*vrough)^2 # Dirichlet energy of vsmooth and vrough"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.4 Convolution\n",
    "The Julia package `DSP` (Digital Signal Processing) includes a convolution function `conv`. After adding this package, the command `conv(a,b)` can be used to compute the convolution of the vectors `a` and `b`. Let’s use this to find the coefficients of the polynomial \n",
    "$$\n",
    "p(x) = (1 + x)(2 − x+ x^2)(1 + x− 2x^2) = 2 + 3x− 3x^2 − x^3 + x^4 − 2x^5.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m  Updating\u001b[22m\u001b[39m registry at `~/.julia/registries/General`\n",
      "\u001b[32m\u001b[1m  Updating\u001b[22m\u001b[39m git-repo `https://github.com/JuliaRegistries/General.git`\n",
      "\u001b[2K\u001b[?25h[1mFetching:\u001b[22m\u001b[39m [========================================>]  99.9 %0.0 %\u001b[36m\u001b[1mFetching:\u001b[22m\u001b[39m [=====>                                   ]  11.5 %                          ]  34.5 %9 % [====================>                    ]  49.7 %\u001b[36m\u001b[1mFetching:\u001b[22m\u001b[39m [=======================>                 ]  57.3 %>            ]  69.7 %4 %>      ]  84.9 %   ]  90.6 %\u001b[36m\u001b[1mFetching:\u001b[22m\u001b[39m [=======================================> ]  96.3 %\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m PDMats ───────────────── v0.9.10\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m OffsetArrays ─────────── v0.11.3\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ShiftedArrays ────────── v1.0.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m MacroTools ───────────── v0.5.3\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m NNlib ────────────────── v0.6.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m StatsModels ──────────── v0.6.7\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m StaticArrays ─────────── v0.12.1\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DSP ──────────────────── v0.5.2\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m QuadGK ───────────────── v2.3.1\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DataAPI ──────────────── v1.1.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m TimeZones ────────────── v0.10.3\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m BinaryProvider ───────── v0.5.8\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Zygote ───────────────── v0.4.1\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m PlotUtils ────────────── v0.6.1\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m EzXML ────────────────── v0.9.5\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m PlotlyJS ─────────────── v0.13.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m PlotlyBase ───────────── v0.3.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DiffEqDiffTools ──────── v1.6.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m IterTools ────────────── v1.3.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m FFMPEG ───────────────── v0.2.4\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Parsers ──────────────── v0.3.10\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Compat ───────────────── v2.2.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DecisionTree ─────────── v0.10.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m TableTraitsUtils ─────── v1.0.1\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m FillArrays ───────────── v0.8.2\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m IterableTables ───────── v1.0.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m RData ────────────────── v0.6.3\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Missings ─────────────── v0.4.3\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m NearestNeighbors ─────── v0.4.4\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Rmath ────────────────── v0.6.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ElasticArrays ────────── v1.0.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Interpolations ───────── v0.12.5\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DocStringExtensions ──── v0.8.1\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Optim ────────────────── v0.19.7\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m NaNMath ──────────────── v0.3.3\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m IRTools ──────────────── v0.3.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Polynomials ──────────── v0.6.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ZygoteRules ──────────── v0.2.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Arpack ───────────────── v0.3.2\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Measures ─────────────── v0.3.1\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m NLSolversBase ────────── v7.5.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m VersionParsing ───────── v1.2.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ArrayInterface ───────── v2.1.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Calculus ─────────────── v0.5.1\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m GaussianProcesses ────── v0.11.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ProgressMeter ────────── v1.2.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m TimerOutputs ─────────── v0.5.3\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m PositiveFactorizations ─ v0.2.3\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m ForwardDiff ──────────── v0.10.7\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m FFTW ─────────────────── v1.0.1\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m FastGaussQuadrature ──── v0.4.1\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m DataStructures ───────── v0.17.6\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m HDF5 ─────────────────── v0.12.5\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m HTTP ─────────────────── v0.8.8\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m GR ───────────────────── v0.44.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m PlotThemes ───────────── v1.0.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m RDatasets ────────────── v0.6.6\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m FileIO ───────────────── v1.2.0\n",
      "\u001b[32m\u001b[1m Installed\u001b[22m\u001b[39m Mustache ─────────────── v1.0.0\n",
      "\u001b[32m\u001b[1m  Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.1/Project.toml`\n",
      " \u001b[90m [717857b8]\u001b[39m\u001b[92m + DSP v0.5.2\u001b[39m\n",
      "\u001b[32m\u001b[1m  Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.1/Manifest.toml`\n",
      " \u001b[90m [7d9fca2a]\u001b[39m\u001b[93m ↑ Arpack v0.3.1 ⇒ v0.3.2\u001b[39m\n",
      " \u001b[90m [4fba245c]\u001b[39m\u001b[93m ↑ ArrayInterface v1.2.1 ⇒ v2.1.0\u001b[39m\n",
      " \u001b[90m [b99e7846]\u001b[39m\u001b[93m ↑ BinaryProvider v0.5.6 ⇒ v0.5.8\u001b[39m\n",
      " \u001b[90m [00ebfdb7]\u001b[39m\u001b[91m - CSTParser v0.6.2\u001b[39m\n",
      " \u001b[90m [49dc2e85]\u001b[39m\u001b[93m ↑ Calculus v0.5.0 ⇒ v0.5.1\u001b[39m\n",
      " \u001b[90m [34da2185]\u001b[39m\u001b[93m ↑ Compat v2.1.0 ⇒ v2.2.0\u001b[39m\n",
      " \u001b[90m [717857b8]\u001b[39m\u001b[92m + DSP v0.5.2\u001b[39m\n",
      " \u001b[90m [9a962f9c]\u001b[39m\u001b[93m ↑ DataAPI v1.0.1 ⇒ v1.1.0\u001b[39m\n",
      " \u001b[90m [864edb3b]\u001b[39m\u001b[93m ↑ DataStructures v0.17.0 ⇒ v0.17.6\u001b[39m\n",
      " \u001b[90m [7806a523]\u001b[39m\u001b[93m ↑ DecisionTree v0.8.3 ⇒ v0.10.0\u001b[39m\n",
      " \u001b[90m [01453d9d]\u001b[39m\u001b[93m ↑ DiffEqDiffTools v1.3.0 ⇒ v1.6.0\u001b[39m\n",
      " \u001b[90m [ffbed154]\u001b[39m\u001b[93m ↑ DocStringExtensions v0.8.0 ⇒ v0.8.1\u001b[39m\n",
      " \u001b[90m [fdbdab4c]\u001b[39m\u001b[93m ↑ ElasticArrays v0.4.0 ⇒ v1.0.0\u001b[39m\n",
      " \u001b[90m [8f5d6c58]\u001b[39m\u001b[93m ↑ EzXML v0.9.4 ⇒ v0.9.5\u001b[39m\n",
      " \u001b[90m [c87230d0]\u001b[39m\u001b[93m ↑ FFMPEG v0.2.2 ⇒ v0.2.4\u001b[39m\n",
      " \u001b[90m [7a1cc6ca]\u001b[39m\u001b[93m ↑ FFTW v0.3.0 ⇒ v1.0.1\u001b[39m\n",
      " \u001b[90m [442a2c76]\u001b[39m\u001b[93m ↑ FastGaussQuadrature v0.4.0 ⇒ v0.4.1\u001b[39m\n",
      " \u001b[90m [5789e2e9]\u001b[39m\u001b[93m ↑ FileIO v1.0.7 ⇒ v1.2.0\u001b[39m\n",
      " \u001b[90m [1a297f60]\u001b[39m\u001b[93m ↑ FillArrays v0.7.0 ⇒ v0.8.2\u001b[39m\n",
      " \u001b[90m [f6369f11]\u001b[39m\u001b[93m ↑ ForwardDiff v0.10.3 ⇒ v0.10.7\u001b[39m\n",
      " \u001b[90m [28b8d3ca]\u001b[39m\u001b[93m ↑ GR v0.41.0 ⇒ v0.44.0\u001b[39m\n",
      " \u001b[90m [891a1506]\u001b[39m\u001b[93m ↑ GaussianProcesses v0.9.0 ⇒ v0.11.0\u001b[39m\n",
      " \u001b[90m [f67ccb44]\u001b[39m\u001b[93m ↑ HDF5 v0.12.3 ⇒ v0.12.5\u001b[39m\n",
      " \u001b[90m [cd3eb016]\u001b[39m\u001b[93m ↑ HTTP v0.8.6 ⇒ v0.8.8\u001b[39m\n",
      " \u001b[90m [7869d1d1]\u001b[39m\u001b[92m + IRTools v0.3.0\u001b[39m\n",
      " \u001b[90m [a98d9a8b]\u001b[39m\u001b[93m ↑ Interpolations v0.12.2 ⇒ v0.12.5\u001b[39m\n",
      " \u001b[90m [c8e1da08]\u001b[39m\u001b[93m ↑ IterTools v1.2.0 ⇒ v1.3.0\u001b[39m\n",
      " \u001b[90m [1c8ee90f]\u001b[39m\u001b[93m ↑ IterableTables v0.11.0 ⇒ v1.0.0\u001b[39m\n",
      " \u001b[90m [1914dd2f]\u001b[39m\u001b[93m ↑ MacroTools v0.5.1 ⇒ v0.5.3\u001b[39m\n",
      " \u001b[90m [442fdcdd]\u001b[39m\u001b[93m ↑ Measures v0.3.0 ⇒ v0.3.1\u001b[39m\n",
      " \u001b[90m [e1d29d7a]\u001b[39m\u001b[93m ↑ Missings v0.4.2 ⇒ v0.4.3\u001b[39m\n",
      " \u001b[90m [ffc61752]\u001b[39m\u001b[93m ↑ Mustache v0.5.13 ⇒ v1.0.0\u001b[39m\n",
      " \u001b[90m [d41bc354]\u001b[39m\u001b[93m ↑ NLSolversBase v7.4.1 ⇒ v7.5.0\u001b[39m\n",
      " \u001b[90m [872c559c]\u001b[39m\u001b[92m + NNlib v0.6.0\u001b[39m\n",
      " \u001b[90m [77ba4419]\u001b[39m\u001b[93m ↑ NaNMath v0.3.2 ⇒ v0.3.3\u001b[39m\n",
      " \u001b[90m [b8a86587]\u001b[39m\u001b[93m ↑ NearestNeighbors v0.4.3 ⇒ v0.4.4\u001b[39m\n",
      " \u001b[90m [6fe1bfb0]\u001b[39m\u001b[93m ↑ OffsetArrays v0.11.1 ⇒ v0.11.3\u001b[39m\n",
      " \u001b[90m [429524aa]\u001b[39m\u001b[93m ↑ Optim v0.19.3 ⇒ v0.19.7\u001b[39m\n",
      " \u001b[90m [90014a1f]\u001b[39m\u001b[93m ↑ PDMats v0.9.9 ⇒ v0.9.10\u001b[39m\n",
      " \u001b[90m [69de0a69]\u001b[39m\u001b[93m ↑ Parsers v0.3.7 ⇒ v0.3.10\u001b[39m\n",
      " \u001b[90m [ccf2f8ad]\u001b[39m\u001b[93m ↑ PlotThemes v0.3.0 ⇒ v1.0.0\u001b[39m\n",
      " \u001b[90m [995b91a9]\u001b[39m\u001b[93m ↑ PlotUtils v0.5.8 ⇒ v0.6.1\u001b[39m\n",
      " \u001b[90m [a03496cd]\u001b[39m\u001b[93m ↑ PlotlyBase v0.2.6 ⇒ v0.3.0\u001b[39m\n",
      " \u001b[90m [f0f68f2c]\u001b[39m\u001b[93m ↑ PlotlyJS v0.12.5 ⇒ v0.13.0\u001b[39m\n",
      " \u001b[90m [f27b6e38]\u001b[39m\u001b[93m ↑ Polynomials v0.5.2 ⇒ v0.6.0\u001b[39m\n",
      " \u001b[90m [85a6dd25]\u001b[39m\u001b[93m ↑ PositiveFactorizations v0.2.2 ⇒ v0.2.3\u001b[39m\n",
      " \u001b[90m [92933f4c]\u001b[39m\u001b[92m + ProgressMeter v1.2.0\u001b[39m\n",
      " \u001b[90m [1fd47b50]\u001b[39m\u001b[93m ↑ QuadGK v2.0.3 ⇒ v2.3.1\u001b[39m\n",
      " \u001b[90m [df47a6cb]\u001b[39m\u001b[93m ↑ RData v0.6.2 ⇒ v0.6.3\u001b[39m\n",
      " \u001b[90m [ce6b1742]\u001b[39m\u001b[93m ↑ RDatasets v0.6.1 ⇒ v0.6.6\u001b[39m\n",
      " \u001b[90m [731186ca]\u001b[39m\u001b[91m - RecursiveArrayTools v0.18.6\u001b[39m\n",
      " \u001b[90m [79098fc4]\u001b[39m\u001b[93m ↑ Rmath v0.5.0 ⇒ v0.6.0\u001b[39m\n",
      " \u001b[90m [1277b4bf]\u001b[39m\u001b[93m ↑ ShiftedArrays v0.5.0 ⇒ v1.0.0\u001b[39m\n",
      " \u001b[90m [90137ffa]\u001b[39m\u001b[93m ↑ StaticArrays v0.11.0 ⇒ v0.12.1\u001b[39m\n",
      " \u001b[90m [3eaba693]\u001b[39m\u001b[93m ↑ StatsModels v0.6.3 ⇒ v0.6.7\u001b[39m\n",
      " \u001b[90m [382cd787]\u001b[39m\u001b[93m ↑ TableTraitsUtils v1.0.0 ⇒ v1.0.1\u001b[39m\n",
      " \u001b[90m [f269a46b]\u001b[39m\u001b[93m ↑ TimeZones v0.10.0 ⇒ v0.10.3\u001b[39m\n",
      " \u001b[90m [a759f4b9]\u001b[39m\u001b[92m + TimerOutputs v0.5.3\u001b[39m\n",
      " \u001b[90m [0796e94c]\u001b[39m\u001b[91m - Tokenize v0.5.6\u001b[39m\n",
      " \u001b[90m [81def892]\u001b[39m\u001b[93m ↑ VersionParsing v1.1.3 ⇒ v1.2.0\u001b[39m\n",
      " \u001b[90m [e88e6eb3]\u001b[39m\u001b[92m + Zygote v0.4.1\u001b[39m\n",
      " \u001b[90m [700de1a5]\u001b[39m\u001b[92m + ZygoteRules v0.2.0\u001b[39m\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m  Building\u001b[22m\u001b[39m EzXML ────→ `~/.julia/packages/EzXML/QtGgF/deps/build.log`\n",
      "\u001b[32m\u001b[1m  Building\u001b[22m\u001b[39m TimeZones → `~/.julia/packages/TimeZones/pjvlM/deps/build.log`\n",
      "\u001b[32m\u001b[1m  Building\u001b[22m\u001b[39m PlotlyJS ─→ `~/.julia/packages/PlotlyJS/b9Efu/deps/build.log`\n",
      "\u001b[32m\u001b[1m  Building\u001b[22m\u001b[39m FFMPEG ───→ `~/.julia/packages/FFMPEG/guN1x/deps/build.log`\n",
      "\u001b[32m\u001b[1m  Building\u001b[22m\u001b[39m Rmath ────→ `~/.julia/packages/Rmath/BoBag/deps/build.log`\n",
      "\u001b[32m\u001b[1m  Building\u001b[22m\u001b[39m Arpack ───→ `~/.julia/packages/Arpack/zCmTA/deps/build.log`\n",
      "\u001b[32m\u001b[1m  Building\u001b[22m\u001b[39m GR ───────→ `~/.julia/packages/GR/oiZD3/deps/build.log`\n",
      "\u001b[32m\u001b[1m  Building\u001b[22m\u001b[39m FFTW ─────→ `~/.julia/packages/FFTW/MJ7kl/deps/build.log`\n",
      "\u001b[32m\u001b[1m  Building\u001b[22m\u001b[39m HDF5 ─────→ `~/.julia/packages/HDF5/Zh9on/deps/build.log`\n"
     ]
    }
   ],
   "source": [
    "# import Pkg; Pkg.add(\"DSP\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Precompiling DSP [717857b8-e6f2-59f4-9121-6e50c889abd2]\n",
      "└ @ Base loading.jl:1186\n",
      "┌ Warning: Module Compat with build ID 120767831084646 is missing from the cache.\n",
      "│ This may mean Compat [34da2185-b29b-5c13-b0c7-acf172513d20] does not support precompilation but is imported by a module that does.\n",
      "└ @ Base loading.jl:947\n",
      "┌ Info: Recompiling stale cache file /Users/vb/.julia/compiled/v1.1/Polynomials/OaK78.ji for Polynomials [f27b6e38-b328-58d1-80ce-0feddd5e7a45]\n",
      "└ @ Base loading.jl:1184\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "6-element Array{Int64,1}:\n",
       "  2\n",
       "  3\n",
       " -3\n",
       " -1\n",
       "  1\n",
       " -2"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using DSP\n",
    "a = [1,1]; # coefficients of 1+x\n",
    "b = [2,-1,1]; # coefficients of 2-x+x^2\n",
    "c = [1,1,-2]; # coefficients of 1+x-2x^2\n",
    "d = conv(conv(a,b),c) # coefficients of product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s write a function that creates a $Toeplitz$ matrix, and check it against the\n",
    "`conv` function. We will also check that Julia is using the very efficient method for\n",
    "computing the convolution.\n",
    "\n",
    "To construct the $Toeplitz$ matrix $T(b)$ defined in equation ([7.3](https://web.stanford.edu/~boyd/vmls/vmls.pdf#equation.7.4.3)) of VMLS, we\n",
    "first create a zero matrix of the correct dimensions $((n + m − 1) × n)$ and then\n",
    "add the coefficients bi one by one. Single-index indexing comes in handy for this\n",
    "purpose. The single-index indexes of the elements $b_i$ in the matrix $T (b)$ are $i,\n",
    "i+m+ n, i+ 2(m+ n), . . . , i+ (n− 1)(m+ n)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6×4 Array{Float64,2}:\n",
       " -1.0   0.0   0.0   0.0\n",
       "  2.0  -1.0   0.0   0.0\n",
       "  3.0   2.0  -1.0   0.0\n",
       "  0.0   3.0   2.0  -1.0\n",
       "  0.0   0.0   3.0   2.0\n",
       "  0.0   0.0   0.0   3.0"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function toeplitz(b,n)\n",
    "    m = length(b)\n",
    "    T = zeros(n+m-1,n)\n",
    "    for i=1:m\n",
    "        T[i : n+m : end] .= b[i]\n",
    "    end\n",
    "    return T\n",
    "end\n",
    "b = [-1,2,3]; a = [-2,3,-1,1];\n",
    "Tb = toeplitz(b, length(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([2.0, -7.0, 1.0, 6.0, -1.0, 3.0], [2, -7, 1, 6, -1, 3])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Tb*a, conv(b,a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.170486 seconds (69.42 k allocations: 64.564 MiB, 40.14% gc time)\n",
      "  0.013096 seconds (165 allocations: 259.266 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.4927701549377119e-12"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = 2000; n = 2000;\n",
    "b = randn(n); a=randn(m);\n",
    "@time ctoep = toeplitz(b,n)*a;\n",
    "@time cconv = conv(a,b);\n",
    "norm(ctoep - cconv)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.1.1",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
